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Chapter 1 
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We review exact approaches and recent results related to the relaxation 
dynamics and description after relaxation of various one-dimensional 
lattice systems of hard-core bosons after a sudden quench. We first 
analyze the integrable case, where the combination of analytical insights 
and computational techniques enable one to study large systems sizes. 
Thermalization does not occur in that regime. However, after relaxation, 
observables can be described by a generalization of the Gibbs ensemble. 
We then utilize full exact diagonalization to study what happens as 
integrability is broken. We show that thermalization does occur in finite 
nonintegrable systems provided they are sufficiently far away from the 
integrable point. We argue that the onset of thermalization can be 
understood in terms of the Eigenstate Thermalization Hypothesis. 

1.1. Introduction 

Understanding how statistical properties emerge from microscopic models 
of many-particle systems has always been of fundamental interest in several 
fields in physics. This topic has been extensively studied in the context of 
classical systems. We know that if we perturb a generic isolated gas in many 
different ways, it will still relax to a unique (Maxwell) velocity distribution 
determined by its energy. This universal behavior (thermalization) has been 
understood in terms of dynamical chaos, namely, the nonlinear equations 
that drive the dynamics ensure that the system explores ergodically all the 
available phase space pQ. However, there is a class of models, known as 
integrable models, for which the presence of a complete set of conserved 
quantities precludes thermalization. In this case, dictated by the initial 
conditions, the dynamics is restricted to a limited region of phase space. 
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More than fifty years ago, Fermi, Pasta, and Ulam (FPU) [5] set up one of 
the first numerical experiments to study how thermalization takes place in 
a one-dimensional (ID) lattice of harmonic oscillators once nonlinear cou- 
plings were added. No signals of ergodicity were found. Those unexpected 
results led to intensive research [3] and ultimately to the development of 
the modern chaos theory [4]. 

Recent advances in cooling and trapping atomic gases has led to in- 
creased interest in understanding what happens in the quantum case. In 
those experiments, the high degree of isolation, combined with the possi- 
bility of controlling interactions and the effective dimensionality of the gas, 
has allowed experimentalists to realize [SHI] and explore the dynamics [HI [H] 
of nearly integrable ID systems. Thermalization was not observed in one 
of the experiments [5] but was indirectly confirmed in the other [S]. Those 
results have motivated intense theoretical research on the dynamics and 
thermalization of isolated quantum systems after a sudden quench, both in 
the integrable [I0H22] and nonintegrable [23H36] regimes. 

Here, we review results for ID systems of hard-core bosons (HCBs) on 
a lattice. We show that thermalization does not occur (in general) when 
the system is integrable. However, observables after relaxation can be de- 
scribed by a generalization of the Gibbs ensemble [TQl [12] . As integrability 
is broken, thermalization does take place [30] EE [36] , and is shown to follow 
after the eigenstate thermalization hypothesis [25j [37l [38] . 

1.2. Methodology 



The HCB Hamiltonian of interest reads 




where t (t r ) is the nearest (next- nearest) neighbor hopping, V (V) is the 
nearest (next-nearest) neighbor interaction, and l^ cxt is an external poten- 
tial. The HCB creation (annihilation) operator in each site is denoted by b\ 
(6J, and the site number occupation by h i = b\b i . We note that in addition 
to the standard commutation relations for bosons, HCB operators satisfy 
the additional constraints b\ 2 = b\ = 0, which preclude multiple occupancy 
of the lattice sites. For t' = V = 0, and any value of V, this model is 
integrable. The approaches used to study this model are described below. 
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1.2.1. Integrable case with V — t' — V — 

This problem can be exactly solved if one realizes that the HCB Hamil- 
tonian can be mapped onto a spin Hamiltonian by means of the Holstein- 
Primakoff transformation |39j . 

4 = S 5 V 1 " 8 ?**' °5 \ 1 '•>, V o? = *ft-l/2, (1.2) 
and that the spin Hamiltonian can be mapped onto a noninteracting fermion 
Hamiltonian utilizing the Jordan- Wigner transformation [40l HI] 

4 = /] ftf* , o-r = e™ ^*<i #A CT | = /t/. _ 1/2. (1.3) 

For simplicity, we will assume open boundary conditions. The resulting 
Hamiltonian for the noninteracting fermions reads 

L-1 L 

H F = -tJ2 (fifi+i + H.c.) + ]T K ex 7//i, (1.4) 
i=i i=i 

and, since it is quadratic, it can be easily diagonalized. Hence, HCBs 
and noninteracting fermions share the same spectrum. The density profiles 
and any density-density correlations will also coincide in both systems. 
The nontrivial differences between HCBs and noninteracting fermions are 
revealed by the off-diagonal correlations. In particular, we will be interested 
in the time evolution of the equal-time one-particle correlations pij , needed 
to compute the momentum distribution function. Once again, using (|1.2[) 

= <SjS_ 7 > = (4°7) = Wj^t + < ^( 1 - 2<J i~4))> anci wc focus 011 now 

to compute Gij = (a^a^). Using (ll.3[) . Gij{r) can be written as [4*21 14*5] 

4-1 3-1 

G«(r) = (* F (r)\ [] j'M'fifl II ^ i7r/!+/! l*F(r)), (1.5) 
fe=i 2=1 

where |^(r)) = e^**/*]**), \9' F ) = UZ^ELi PnJl\0) is the 
initial state (a Slater determinant), N is the number of particles, and 
r the time variable. The action of exponentials whose exponents arc 
bilinear on fermionic creation and annihilation operators, as Hp is in 
this case, on Slater determinants generate new Slater determinants, so 
I*f(t)) = n^=iELi P nm(r)/t|0). The matrix P(r) can be computed 
as P(t) = e-^^P 7 = Ue- <TE / fi Utp J , where we have used that 
Hi?U = UE, where E is a diagonal matrix containing the eigenenergies 
and U is the unitary matrix of eigenvectors. Furthermore, the action 
of n?=i eT w ^ 1 ' 1 on \^f(t)) changes the sign of P nm (r) for n < j — 1, 
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m = 1, . . . , N, and the creation of a particle at site j implies the addi- 
tion of a column with only one nonzero element (Pj n+i(t) = 1) [the same 
applies to the action of nl=i e ln ^fk fa on the left of Eq. (|1.5|) ]. Hence, 

N L N L 

Ga{r) = <0| J! E Ktn(r)f n II E 1°)' (!-6) 

m=l n=l fc=l J=l 

= det[(P < (r)) t I" , (r)]. (1.7) 
In Eq. (|1.6|) . the matrix elements P^ m (r) and Pi k (r) have the form 

r -P/3 7 (r) for /3< a, 7 = !,•••, # 
i^ 7 (r) = { P^(r)for^>a, 7 = l,...,iV , (1.8) 
[ S a p for 7 = N + 1 



with a = /3 = n, I, and 7 = m, fc. Equation (|1.7p follows from (|1.6[) by 
using the identity 

<0|/ ai • • • f aN Jl N+l ■ ■ ■ fl\0) = e Xl - XN+1 S ai p Xl ■ ■ ■ S aN+1 p XN+i , (1.9) 

where £ a i - a n+i j s ^ ne Levi-Civita symbol in N + 1 dimensions, and the 
indices A have values between 1 and N + 1. Employing Eq. (|1.7[) . can 
be calculated in polynomial time, which scales as L 2 N 3 , using a computer. 

We will also be interested in describing the momentum distribution 
function after relaxation by using statistical ensembles. A polynomial time 
approach in this case is only known within the grand-canonical formalism 
[44] . The one-particle density matrix in this ensemble can be written as 

p l3 = -Tr \b\b 3 e -bt j , Z = Tr. 

pa = ^Tr n em4/fce ^ n e_M/! /! . 
1 k=i 1=1 ) 

where, /i is the chemical potential, fee is the Boltzmann constant, T is 
the temperature, and Z (which is the same for HCBs and fermions) is 
the partition function. To arrive to Eq. (|1 . 10|) . in addition to the Jordan- 
Wigner transformation (| 1 . 3[) , we have used the cyclic property of the trace. 
Another property of the trace, this time over the fermionic Fock space, is 

(1.11) 

where I is the identity matrix. Equation (jl.lll) allows us to compute Z 
as Z = Y\i [l + e _ ^" _A1 ^ fcsT ] . By noticing that for i =/= j, we can write 



x y zn 

e e • • • e 
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fjfj = exp (J2 mn fL A mnfn) ~ 1) where the only nonzero element of A is 
1, the off-diagonal elements of pij (i ^ j) can be obtained as 
1 



Aij 



Pij 



Z 

-det 



^det 1 1 + (I + A)OiUe- (E - AlI)/A;BT U 1 '02 
I + 1 Ue-( E -^/ fe - T Ut0 2 ]}, 



(1.12) 

where Oi (O2) is diagonal with the first j — 1 (i— 1) elements of the diagonal 
equal to —1 and the others equal to 1. The diagonal elements of pij are the 
same as for noninteracting fermions and can be computed as 

-1 



Pii = 



,-(H F - M I)/fc B T 



1 _1 








a 





,-(B-rt/i B T 



ut 



(1.13) 



The computational time within this approach scales as L 5 . 

The momentum distribution function n(k), in and out of equilibrium, 
is then determined by the expression n(k) = J^i 



-ifc(Z-m) 
. e yral • 



1.2.2. Nointegrable case with V 



ext 







For this case, we will make use of full exact diagonalization (see e.g., [45]). 
This approach has the disadvantage that the dimension of the matrices 
that need to be diagonalized scales exponentially with system size. Since 
the Hamiltonian (jl.ljl conserves the total number of particles, we work with 
a fixed number of particles N — L/3, which reduces the dimension of our 
problem from 2 L to (j^) . To further reduce the dimension of the matrices to 
be diagonalized, we consider systems with periodic boundary conditions and 
no external potential (V cxt — 0). Then, by using translational symmetry, 
we can block diagonalize the Hamiltonian, with the size of each momentum 
block being ~ 1/L the size of the original matrix. All momentum sectors, 
whose dimensions are shown in the table below, are diagonalized. They are 
all used to construct the microcanonical and canonical ensembles. 



Dimensions of all momentum sectors 


L = 18 


k = 0,6 


fc = 1,5,7 


fc = 2,4,8 


fc = 3,9 


dimension 


1038 


1026 


1035 


1028 


L = 21 


fc = 0,7 


other fc's 






dimension 


5538 


5537 






L = 24 


fc = 0,8 


fc = 4,12 


fc = 2,6,10 


odd fc's 


dimension 


30667 


30666 


30664 


30624 
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1.3. Results 

We focus on the dynamics after a sudden quench. This means that we 
start with some eigenstate of an initial Hamiltonian, which may not be the 
ground state, then at r = some parameter is changed and the system 
is allowed to evolve. Independently of whether the Hamiltonian is inte- 
grable or not, one can always write the initial state wavef unction IV'ini) 
in the eigenstate basis of the final Hamiltonian, i.e., |V>ini) = ^a^al^a) 
with C Q = (vf'alV'ini) and H\^ a ) = E a \& a ). The dynamics of the wave- 
function takes the form |V>(f)) = e~ iTA/h \^ixA) = T, a e~ iTEa/h C a \^f a ) and 
the expectation value of any observable O can be written as {0(t)) = 
m)\0\m) = Ea/^^e'^-^l'O^, where O a/3 = (* a \0\*p). If 
the spectrum is nondegenerate and incommensurate, the infinite time av- 
erage and the observable after relaxation is determined by 

{0) = O diag = Y J \C a \ 2 O aa . (1.14) 

ot 

This exact result can be thought as the prediction of a "diagonal ensem- 
ble", where \C a \ 2 is the weight of each state [23], and is different from any 
conventional ensemble of statistical mechanics. 

1.3.1. Integrable case with V — t' — V — 

Here, our set up is close in spirit to the one in the experiment [8 . The initial 
state is the ground state of a harmonic trap with a staggered potential and, 
at t = 0, we turn off the staggered potential and allow the system to evolve 
in the presence of the trap [12] , In addition to density profiles and n(k), we 
also study the occupation of the natural orbitals, which are the eigenstates 
of the one-particle density matrix, determined by the eigenvalue equation 
Pij4>] — \4>7- The lowest natural orbital is the highest occupied one. 
Figure [TTT1 depicts the evolution of the occupation of the zero-momentum 
state n(k = 0) and the lowest natural orbital Ao when, (i) the initial state 
has a half-filled insulator in the center of the trap [Fig. I1.2f a)] and, (ii) 
two insulating shoulders surround a central superfluid region [Fig. I1.2l d)]. 
In both cases, the two observables undergo relaxation dynamics, which 
ultimately brings them to an almost time independent result. This shows 
that relaxation is not precluded by integrability, and the question that 
remains to be answered is how to describe those observables after relaxation. 
As seen in Fig. 11.11 they are clearly different from the predictions of the 
grand-canonical ensemble (GE in the figures), which are obtained after 
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2500 5000 7500 10000 
X 



2500 5000 7500 1000 
X 



Fig. 1.1. Dynamics of n(k = 0) (top plots) and An (bottom plots) after a staggered 
potential is turned off in harmonically confined systems with 900 lattice sites and a trap 
curvature V% = 3 X 10~ 5 t. r is given in units of h/t, and the evolution starts from the 
ground state in the presence of a staggered potential of strength 0.5t. The number of 
particles is: (a) N = 200 and (b) N = 299. [The corresponding initial density profiles 
can be seen in Figs. [L2f a) and ll.2f dH. The dashed-doted lines depict the results within 
the grand-canonical ensemble (GE), with (a) T = 0.31t and (b) T = 0.334, and the 
dashed lines the results within the generalized Gibbs ensemble (GGE). 



determining the temperature and chemical potential so that 
E = IlY {# e -<*-M*)/*BTj , N = ±Tr{Ne-^-^/ k * T }, (1.15) 

where N = J^. b\b i , and E and N are the energy and number of particles 
in the time evolving state, which are conserved during the evolution. We 
note that for the system sizes considered, finite size effects are negligible. 

The lack of relaxation to the thermal state may not be surprising consid- 
ering that the system is integrable, and, hence, the existence of conserved 
quantities may preclude thermalization. In Ref. |10j , a generalization of the 
Gibbs ensemble (GGE) was proposed in order to account for the conserved 
quantities and still be able to statistically describe integrable systems. The 
density matrix for the GGE was determined by maximizing the many-body 
Gibbs entropy S = fcsTr [p c ln(l/p c )] subject to the constraints imposed by 
all the integrals of motion. The result reads 

/ 3 c = -l e -E i A J / s £ c = T r{e-£' A ^} (1.16) 

where Z c is the generalized partition function, {/;} is a full set of integrals 
of motion, {A;} are the Lagrange multipliers, and I = 1, . . . ,L. The La- 
grange multipliers are computed using the expectation values of the full set 
of integrals of motion in the initial state so that (Ii)im = Tr{/ ZJ o c }. For 
HCBs, which can be mapped to noninteracting fermions, a complete set 
of integrals of motion is provided by the projection operators to the non- 
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* ka r| 

Fig. 1.2. Initial state and time average values of: (a),(d) density profiles, (b),(e) mo- 
mentum distribution functions, and (c),(f) occupation of the lowest 100 natural orbitals. 
The averages are computed between r = 5000/i/t and r = 10000/l/i with measurements 
done in time intervals At = 40i, and correspond to the dynamics depicted in Fig. 11.11 
The results of the time average are compared with those obtained in the grand-canonical 
ensemble (GE) and the generalized Gibbs ensemble (GGE) explained in the text. The 
number of particles is N = 200 (a)-(c) and N = 299 (d)-(f). In (a) and (d), for the 
initial state, the occupations plotted are the averaged density per unit cell. Note that in 
the presence of the staggered potential, the density exhibits large fluctuations from site 
to site. Flat regions of the unit cell occupations correspond to insulating domains 1121 . 




interacting single particle eigenstates {/;} = {7/7/}, where {7; } ({7/}) 
creates (annihilates) a single particle in an eigenstate of Eq. (| 1 .4[) . The 
resulting Lagrange multipliers read A; = ln[(l — {^/}im)/(^)ini]- They allow 
one to build the density matrix in Eq. (j!.16[) and to compute expectation 
values as it was described for the grand-canonical ensemble in Sec. 11.2.11 

Figure im shows that the GGE calculations for n(k = 0) and Ao properly 
predict the outcome of the relaxation dynamics. We have also computed 
the time average (between r = 5000fi/i and lOOOOh/t) of the full density 
profiles, n(k), and A,,. They are shown in Fig. 11.21 There, the time averages 
are compared with the results for the initial state and with the predictions 
of the GE and the GGE. That comparison clearly shows that, unlike the 
GE, the GGE is able to predict all those single particle observables after 
relaxation. Note that, when written in the bosonic language, the constraints 
lose the bilinear character they have in the fcrmionic representation, i.e., the 
outcome of the GGE calculation is not at all trivial, as it would be if done 
for noninter acting fermions. Recent numerical and analytical studies have 
addressed various aspects of the GGE pDHTOl ITS] . However, why the GGE 
and the diagonal ensemble should agree still remains to be understood. 
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1.3.2. Nonintegrable case with V 



ext 







To study the effects of breaking integrability, we prepare an initial state that 
is an eigenstate of a Hamiltonian (in the total momentum k = sector) with 
tini, Vi n i, t', V and then quench the nearest-neighbor parameters to tj in , 
Vfi n without changing t', V' . The same quench is repeated for different 
values of t' , V as one departs from t' = V = [30] . To find whether the 
dynamics brings the observables to the predictions of the diagonal ensemble 
(|1.14p . we calculate the normalized area between the observables during the 
time evolution and their infinite time average, i.e., 5n^(r) — (J\ \ n (k, T ) ~ 
n diag(k)\)/Y^,k ndiag{k). Similarly, we compute 5Nk for the structure factor 
N(k), which is the Fourier transform of the density-density correlations. 

In Fig. 11.31 we show results for 5rik and 5Nk vs r for three different 
quenches and two system sizes. The time evolution is very similar in all 
cases, and is consistent with a fast relaxation of both observables towards 
the diagonal ensemble prediction (in a time scale r ~ h/t). The average 
differences after relaxation and their fluctuations can be seen to decrease 
with increasing system size. From these results, we infer that, for very 
large systems sizes, n(k) and N(k) should in general relax to exactly the 
predictions of Eq. (II . 14)) even if the system is very close or at integrability. 

We then say that thermalization takes place if the results of conventional 
statistical ensembles and those of the diagonal ensemble are the same. In 
Fig. 11.4( a). we compare the diagonal ensemble results with the predictions 



o.i 







t-=y=0.24 



— ' 24 



t'=V'=0.03 





. 0.1 



t'=V=0 



■ L=21 

■ L=24 



I - ^ ' — i — i — i — i — i 



0.1 



to 



0.1 





0.1 



20 40 60 80 100 20 40 60 80 100 
X X 



Fig. 1.3. Evolution of Sn^ (left panels) and 8Nk (right panels) after a quench from 
t ini = 0.5t, V ini = 2.0t tot/in = t, Vfi n = t, with t' ini = t' fin = t' and V[ nl = V} m = V , 
for two system sizes. The initial state was selected within the eigenstates with total 
momentum k = such that after the quench the effective temperature is T = 3.0 in all 
cases. Given the energy of the initial state E, T follows from E = Z- 1 Tr{He- i */ fc s T }, 
where Z = Tr{e — H / k B T }. The trace runs over the full spectrum. 
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0.125 0.25 
t'.V 



t'=V=0.24 


i 


i 1 i 
r*=y«=o.a? ■ 


- (b) 2 












i-10 


10 



-15 



-10 



10 



Fig. 1.4. (a) Differences between the predictions of the diagonal and microcanonical 
ensembles (calculated as Sn^ and SN^ in Fig. 11.31) . Results are shown for T = 2.0 and 
T = 3.0. (b) n(k = 0) as a function of the energy for all the eigenstates of the Hamiltonian 
(including all momentum sectors), (main panel) t = V = 1 and t' = V' = 0.24. (inset) 
t = V = 1 and t' = V = 0.03. The systems in (a) and (b) have L = 24 and JV 6 = 8. 



of the microcanonical ensemble for our two observables of interest. Far from 
intcgrability the differences are small and decrease with increasing system 
size [30j . i.e., thermalization takes place. As one approaches integrability, 
the differences increase, signaling a breakdown of thermalization in ID. 

Thermalization away from integrability, as well as its failure close to 
integrability, can be understood in terms of the eigenstate thermalization 
hypothesis (ETH) [25j EH [38]. ETH states that, for generic systems, the 
fluctuations of eigenstate expectation values of observables is small between 
eigenstates that are close in energy, which implies that the microcanonical 
average is identical to the prediction of each eigenstate, which is the same 
as saying that the eigenstates already exhibit thermal behavior. If this 
holds, thermalization in an isolated quantum system will follow for any 
distribution of | C a | 2 that is narrow enough in energy. 

The main panel in Fig. I1.4f b) depicts n(k = 0) [similar results were 
obtained for n(k ^ 0) and N(k)} in each eigenstate of the Hamiltonian when 
the system is far from integrability. After a region of low energies where the 
eigenstate expectation values exhibit large fluctuations, one can see another 
region where fluctuations are small and ETH holds. The inset shows that 
for a system close to the integrable point, in which thermalization is absent 
[Fig. |1.4f a)]. the eigenstate to eigenstate fluctuations of n(k = 0) are very 
large over the entire spectrum and ETH does not hold. 
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